{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import numpy as np"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## algorithm"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def simplex(c, A, b):\n",
    "    table = initialize(c, A, b)\n",
    "    while not search_optimum(table):\n",
    "        pass\n",
    "    return solution(c, table)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def initialize(c, A, b):\n",
    "    (m, n), k = A.shape, len(c)\n",
    "\n",
    "    # simplex table:\n",
    "    # |A|E|b|\n",
    "    # |c|0|0|\n",
    "    table = np.zeros((m + 1, m + n + 1))\n",
    "    table[:m, :n] = A\n",
    "    table[range(m), range(n, n + m)] = 1\n",
    "    table[:-1, -1] = b\n",
    "    table[-1, :k] = c\n",
    "\n",
    "    return table"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def search_optimum(table):\n",
    "    index = np.argwhere(table[-1, :-1] > 0).ravel()\n",
    "    \n",
    "    # optimum found\n",
    "    if not len(index):\n",
    "        return True\n",
    "    \n",
    "    # pivotal column\n",
    "    j = index[0]\n",
    "    column = table[:-1, j].copy()\n",
    "    column[column <= 0] = -1\n",
    "    \n",
    "    if np.all(column <= 0):\n",
    "        raise ArithmeticError('the system is unbounded')\n",
    "\n",
    "    # pivotal row\n",
    "    pivots = table[:-1, -1] / column\n",
    "    pivots[column <= 0] = np.inf\n",
    "    i = np.argmin(pivots).ravel()[0]\n",
    "\n",
    "    # eliminate by pivot at (i, j)\n",
    "    row = table[i] / table[i][j]\n",
    "    table[:] -= np.outer(table[:, j], row)\n",
    "    table[i, :] = row\n",
    "    table[:, j] = table[:, j].round()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def solution(c, table):\n",
    "    (m, n), k = table.shape, len(c)\n",
    "\n",
    "    # pivotal columns\n",
    "    s = np.sum(table == 0, axis=0) == m - 1\n",
    "    t = np.sum(table == 1, axis=0) == 1\n",
    "\n",
    "    # solution\n",
    "    x = np.zeros(n - 1)\n",
    "\n",
    "    for j in range(n - 1):\n",
    "        if s[j] and t[j]:\n",
    "            x[j] = table[:, j] @ table[:, -1]\n",
    "\n",
    "    return dict(\n",
    "        x=x[:k],\n",
    "        slack=x[k:],\n",
    "        max=-table[-1, -1],\n",
    "        table=table,\n",
    "    )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## linear program #1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```\n",
    "maximize: -x + 3y + 2z\n",
    "\n",
    "subject to:\n",
    "x + y + z <= 6\n",
    "x     + z <= 4\n",
    "    y + z <= 3\n",
    "x + y     <= 2\n",
    "\n",
    "x, y, z >= 0\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "c = np.array([-1, 3, 2])\n",
    "A = np.array([\n",
    "    [1, 1, 1],\n",
    "    [1, 0, 1],\n",
    "    [0, 1, 1],\n",
    "    [1, 1, 0],\n",
    "])\n",
    "b = np.array([6, 4, 3, 2])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "x \n",
      " [ 0.  2.  1.] \n",
      "\n",
      "slack \n",
      " [ 3.  3.  0.  0.] \n",
      "\n",
      "table \n",
      " [[ 1.  0.  0.  1.  0. -1.  0.  3.]\n",
      " [ 2.  0.  0.  0.  1. -1.  1.  3.]\n",
      " [-1.  0.  1.  0.  0.  1. -1.  1.]\n",
      " [ 1.  1.  0.  0.  0.  0.  1.  2.]\n",
      " [-2.  0.  0.  0.  0. -2. -1. -8.]] \n",
      "\n",
      "max \n",
      " 8.0 \n",
      "\n"
     ]
    }
   ],
   "source": [
    "lp = simplex(c, A, b)\n",
    "\n",
    "for k in ['x', 'slack', 'table', 'max']:\n",
    "    print(k, '\\n', lp[k], '\\n')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## linear program #2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```\n",
    "maximize: 2r + 4s + 3t + u\n",
    "\n",
    "subject to:\n",
    "3r +  s +  t + 4u <= 12\n",
    " r - 3s + 2t + 3u <= 7\n",
    "2r +  s + 3t -  u <= 10\n",
    "\n",
    "r, s, t, u >= 0\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "c = np.array([2, 4, 3, 1])\n",
    "A = np.array([\n",
    "    [3, 1, 1, 4],\n",
    "    [1, -3, 2, 3],\n",
    "    [2, 1, 3, -1]\n",
    "])\n",
    "b = np.array([12, 7, 10])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "x \n",
      " [  0.   10.4   0.    0.4] \n",
      "\n",
      "slack \n",
      " [  0.  37.   0.] \n",
      "\n",
      "table \n",
      " [[  0.2   0.   -0.4   1.    0.2   0.   -0.2   0.4]\n",
      " [  7.    0.   11.    0.    0.    1.    3.   37. ]\n",
      " [  2.2   1.    2.6   0.    0.2   0.    0.8  10.4]\n",
      " [ -7.    0.   -7.    0.   -1.    0.   -3.  -42. ]] \n",
      "\n",
      "max \n",
      " 42.0 \n",
      "\n"
     ]
    }
   ],
   "source": [
    "lp = simplex(c, A, b)\n",
    "\n",
    "for k in ['x', 'slack', 'table', 'max']:\n",
    "    print(k, '\\n', lp[k], '\\n')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
